home *** CD-ROM | disk | FTP | other *** search
/ POINT Software Programming / PPROG1.ISO / pascal / swag / math.swg / 0075_FFT algorithm.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1994-08-25  |  2.3 KB  |  103 lines
  1. {
  2. From: marcel.hoogeveen@hacom.wlink.nl (Marcel Hoogeveen)
  3.  
  4. GR> FFT stands for Fast Fourier Transform.  It is a quick way to conver
  5. GR> time domain data (ie oscilliscopy data with time on the x-axis) to
  6. GR> frequency domain (frequency on the x-axis, like a frequency spectrum
  7. GR> analyzer).  This is a usefull data analysis method.  I would also like
  8. GR> to get some source for this.
  9.  
  10.  
  11. This is what i have of FFT source code, it should work if you tweak it a bit.
  12. (It did for me when i used it in my analasis program).
  13. Don't ask me how it works, i know how a DFT works but a FFT well .. just use
  14. the source. :)
  15.  
  16. }
  17. Program FFT;
  18. Const Twopi=6.283185303;
  19.  
  20. Type Curve=array[1..nfft] of real;
  21.  
  22. Var {This is for you to find out}
  23.  
  24. { Calculation of the Discrete Fourier Transfor }
  25. { Using a Fast Fourier Transform algorithm     }
  26. {                                              }
  27. { XR and XI are array of reals !!!             }
  28. { They contain on entry the input sequence and }
  29. { on return the transfrom                      }
  30. { ISI defines the transform direction          }
  31. { If ISI=-1 then forward, if ISI=1 then invert }
  32. {                                              }
  33. { The dimension is 2**M                        }
  34.  
  35. Procedure RFFT (VAR XR,XI:Curve;  N:integer;  ISI:Integer);
  36. Var
  37. M,NV2,LE,LE1,IP,I,J,K,L: Integer;
  38. C,THETA,UR,UI,TR,TI:Real;
  39.  
  40. Begin
  41. M:=Round(LN(N)/LN(2));
  42. NV2:= N DIV 2;
  43. J:=1;
  44. For I:= 1 to N-1 do
  45. Begin
  46. If (I<J) then
  47. Begin
  48. TR:=XR[J];            TI:=X[J];
  49. XR[J]:=XR[I];         XI[J]:=XI[I];
  50. XR[I]:=TR;            XI[I]:=TI;
  51. End;
  52. K:=NV2;
  53. While (K<J) do
  54. Begin
  55. J:=J-K;
  56. K:=K DIV 2;
  57. End;
  58. J:=J+K;
  59. End;
  60. LE:=1;
  61. C:=ISI*TWOPI;
  62. For L:=1 TO M do
  63. Begin
  64. LE1:=LE;
  65. LE:=LE*2;
  66. For J:=1 TO LE1 do
  67. Begin
  68. THETA:= C*(J-1)/LE;
  69. UR:=COS(tHETA);
  70. UI:=SIN(THETA);
  71. I:=J;
  72. Repeat
  73. IP:=I+LE1;
  74. TR:=XR[IP]*UR-XI[IP]*UI;
  75. TI:+XR[IP]*UI+XI[IP]*UR;
  76. XR[IP]:=XR[I]-TR;           XI[IP]:=XI[I]-TP;
  77. XR[I]:=XR[I]+TR;            XI[I]:=XI[I]+TI;
  78. I:=I+LE;
  79. Until (I>=N)
  80. End;
  81. End;
  82. If ISI=-1 then
  83. Begin
  84. For I:= 1 TO N do
  85. Begin
  86. XR[I]:=4*XR[I]/N;             XI[I]:=4*XI[I]/N;
  87. End;
  88. End;
  89. End;
  90.  
  91.  
  92. Begin
  93. For I := 1 to NUMSAM do
  94. Begin
  95. FREAL[I]:=SAMPLEBUFFER[I];
  96. FIMAG[I]:=0;
  97. End;
  98. RFFT(FREAL,FIMAG,NUMSAM,-1);
  99. DC:=FREAL[1]/2;
  100. For I:= 1 to NUMSAM dO
  101. FREAL[I]:=FREAL[I]*FREAL[I]+fIMAG[I]*FIMAG[I];
  102. End.
  103.